function[raw_data] = getCOM(videoname, varargin)
% Function to perform the initial steps of video analysis: detects
% molecules and finds their COMs, intensities, and Rgs. Input should be the
% name of a .fits file (in the MATLAB directory). Optionally, the frames to
% be analyzed can be specified as well: 
%
%   getCOM('foobar.fits') analyzes the entire movie
%
%   getCOM('foobar.fits', 123) analyzes frame 123 of the movie
%
%   getCOM('foobar.fits', 123, 456) analyzes frames 123 to 456
%
%   getCOM('foobar.fits', 123, 'hello world') analyzes frame 123 and saves 
%   'hello world' as a comment in the .mat file
%   
%   getCOM('foobar.fits', 123, 456, 'hello world') analyzes frames 123 to
%   456 and saves the comment 'hello world'
%   
%   getCOM('foobar.fits', 123, 456, 111, 222, 333, 444) analyzes frames 123
%   to 456 in the region x = 111~222 and y = 333~444.
%   
%   getCOM('foobar.fits', 123, 456, 111, 222, 333, 444, letter) analyzes 
%   frames 123 to 456, in the region x = 111~222 and y = 333~444, and 
%   appends the character 'letter' to the end of the saved .mat file.
%
%   getCOM('foobar.fits', 123, 456, 111, 222, 333, 444, letter, save_dir) 
%   analyzes frames 123 to 456, in the region x = 111~222 and y = 333~444,
%   appends the character 'letter' to the end of the saved .mat file, and
%   saves the .mat files into the directory defined by save_dir. 
%   NOTE 1 SEP 2014: If save_dir is not given, the .mat file will be save 
%   in the present working directory (which, from call_getCOM.m is 
%   currently set BlueWolves.)
%
%   getCOM('foobar.fits', 123, 456, 111, 222, 333, 444, letter, save_dir) 
%   analyzes frames 123 to 456, in the region x = 111~222 and y = 333~444,
%   appends the character 'letter' to the end of the saved .mat file,
%   saves the .mat files into the directory defined by save_dir, and saves 
%   the comment 'hello world'


% Clear all variables except videoname & varargin.
clearvars -except videoname varargin

% Parse getCOM_Master arguments.
[min, max, xmin, xmax, ymin, ymax, letter, save_dir, comments] = ...
                                                parse_inputs(varargin{:});
                                     
% Change directory to save_dir. This is the directory where the .fits
% video resides.
cd(save_dir);

% Read video.
video = fitsread(videoname);
video = double(video);

% Extract the KCT from the header data.
info = fitsinfo(videoname);
header = info.PrimaryData.Keywords;
KCT = header{find(strcmp(header(:,1),'KCT')),2};

% Change videoname to format used for storage (no path, no .fits).
if(~isempty(strfind(videoname, '/')))
    slashLocs = strfind(videoname, '/');
    videoname = videoname(slashLocs(size(slashLocs, 2)) + ...
                                                    1:size(videoname, 2));
end
videoname = videoname(1:size(videoname, 2) - 5);

% Analyze all frames of the video if no min/max frames are given (see
% function parse_inputs for details.)
if(max == -1)
    max = size(video, 3);
end

% Set analysis settings.
% Video name
name = videoname;
% Number of frames
num_frames = size(video, 3);
% Number of frames analyzed
num_frames_analyzed = max - min + 1;
% Number of frames added over for clustering analysis, centered on current
% frame. For example, num_accumulated = 1 means only 1 frame is being used,
% num_accumulated = 3 means one frame on either side. Should always be odd.
num_accumulated = 1;
% Kernel size of Gaussian used in initial smoothing
gaussian_kernel_size = 3;
% Standard deviations of Gaussian used in initial smoothing
gaussian_std = 1.0;
% Numer of RMS' above the mean the threshold is
threshold_RMS_multiplier = 4.0;

% UPDATE 9 DEC 2014: Error check if threshold_RMS_multiplier ~= 5.0.
% if threshold_RMS_multiplier ~= 5.0
%     'THE VARIABLE, threshold_RMS_multiplier, is not 5.0! '
%     'PROGRAM PAUSED! PRESS ANY KEY TO CONTINUE, OR ABORT PROGRAM.'
%     pause
% end
    
% Object used for morphological closing of clustering
close_object = strel('disk', 2);
% Object used for morphological opening of clustering
open_object = strel('disk', 2);
% Object used for dilating clustering
dilate_object = strel('disk', 2);
% Object used for morphological closing of the edges
edge_close_object = strel('disk', 2);
% Object used for morphological opening of the edges
edge_open_object = strel('disk', 2);
% Object used for dilating edge detection
edge_dilate_object = strel('disk', 2);
% Number of data points (for raw_data), must be manually changed!
num_data = 7;
% Any comments (replaced by comments feature in parse_inputs)
% comments = '';

% For Python importing, the strels must be turned into arrays:
close_object_array = getnhood(close_object);
open_object_array = getnhood(open_object);
dilate_object_array = getnhood(dilate_object);
edge_close_object_array = getnhood(edge_close_object);
edge_open_object_array = getnhood(edge_open_object);
edge_dilate_object_array = getnhood(edge_dilate_object);

% Start of program, using helper functions defined below.
for place_in_series = min:max
    
    % Accumulate frames for clustering analysis. Centered on current frame.
    image = accumulate_frames(video, num_frames, place_in_series, num_accumulated); %accumulates 
    
    % Find background intensity. Initial crude approximation: divide total
    % intensity by total number of pixels.
    background = initial_background_intensity(image);
    
    % Subtract the background from the original image.
    image2 = background_subtract(image, background);
    
    % Use a Gaussian filter to smooth the image.
    smoothed = gaussian_filter(image2, gaussian_kernel_size, gaussian_std);
    
    % Find the RMS intensity of the image.
    rms = initial_background_rms(smoothed);
    
    % Detect pixels that are above threshold_RMS_multiplier * RMS value.
    BWdilated{place_in_series} = clustering_detection(smoothed, rms, threshold_RMS_multiplier, close_object, open_object, dilate_object);
    
    % Refinement of background calculation based on first pass of clustering
    % analysis.
    background = calculate_background_intensity(image, BWdilated{place_in_series});
    
    % Subtract the newly calculated background intensity from the ORIGINAL
    % image. Set the lower limit on the intensity value to be zero
    % (pixel = 0).
    image2 = background_subtract(image, background);
    
    % Use a Gaussian filter to smooth the image.
    smoothed = gaussian_filter(image2, gaussian_kernel_size, gaussian_std);
    
    % Find the RMS intensity of the image.
    rms = calculate_background_rms(smoothed, BWdilated{place_in_series});
    
    % Set bounds. Molecule must be inside these pixel values in each frame.
    for y = 1:size(smoothed, 1)
        for x = 1:size(smoothed, 2)
            if(y < ymin || y > ymax)
                smoothed(y, x) = 0;
                image2(y, x) = 0;
            elseif(x < xmin || x > xmax)
                smoothed(y, x) = 0;
                image2(y, x) = 0;
            end
        end
    end
    
    % Repeat clustering algorithms.
    BWdilated{place_in_series} = clustering_detection(smoothed, rms, threshold_RMS_multiplier, close_object, open_object, dilate_object);
    
    % Execute Jaeyoon's edge detection function.
    [edges, edgeMap{place_in_series}, connected, number_of_molecules] = edge_detection(image2, BWdilated{place_in_series}, edge_close_object, edge_open_object, edge_dilate_object);
    
    % Execute Jaeyoon's molecule extension length function.
    extension_lengths = extension_length(edges, connected, number_of_molecules);
    
    % Catalog the data.
    data = catalog(image2, connected, number_of_molecules);
    
    % Imshow the original image and the connected image. Comment out if you
    % don't need to display these figures.
    %imshow(mat2gray(image));
    %pause;
    %imshow(mat2gray(image) + connected*0.2); 
    %pause;

    % Create the variable raw_data. This is what gets sent to my Python
    % script for archiving and further analysis.
    for i = 1:number_of_molecules
        for j = 1:num_data - 1
            raw_data(place_in_series - min + 1, i, j) = data(i,j);
        end
        raw_data(place_in_series - min + 1, i, num_data) = extension_lengths(i);
    end
    if(number_of_molecules == 0)
        raw_data(place_in_series - min + 1, 1, 1:num_data - 1) = zeros(1, num_data - 1);
        raw_data(place_in_series - min + 1, 1, num_data) = 0;
    end
end

% filter "molecules" that are too dim/too bright to consider.
raw_data = filter_data(raw_data);

% Save this .fits video's raw_data xCOM and yCOM to a text file for error 
% checking.
output_name = strcat(name, letter);
QAcheckfile = strcat(name, 'getCOMCheck.txt');
fileID = fopen(QAcheckfile, 'a');
fprintf(fileID, strcat(output_name, ':\n'));
formatSpec = strcat('Start frame is %i, End frame is %i.\nxmin, xmax, ', ...
                    'ymin, ymax are %i, %i, %i, %i.\n');
fprintf(fileID, formatSpec, min, max, xmin, xmax, ymin, ymax);
fprintf(fileID, '\nxCOM \n');
dlmwrite(QAcheckfile, raw_data(:, :, 1), '-append');
fprintf(fileID, '\nyCOM \n');
dlmwrite(QAcheckfile, raw_data(:, :, 2), '-append');
fprintf(fileID, '--------------------------\n');
fclose(fileID);

% Save the .mat files. Set variable 'name' to 'output_name'. This is 
% important for subsequent MasterDataFrame.py processing.
name = output_name;
save(strcat(save_dir, '/', output_name), 'name', 'num_frames' , 'num_accumulated', ...
    'gaussian_kernel_size', 'gaussian_std', 'threshold_RMS_multiplier', ...
    'close_object_array', 'open_object_array', 'dilate_object_array', ...
    'edge_close_object_array', 'edge_open_object_array', ...
    'edge_dilate_object_array', 'num_data', 'comments', 'raw_data', 'KCT');

BWdilated = cell2mat(BWdilated);
BWdilated = reshape(BWdilated, 512, 512, max-min+1);
BWdilated = permute(BWdilated, [3 1 2]);
edgeMap = cell2mat(edgeMap);
edgeMap = reshape(edgeMap, 512, 512, max-min+1);
edgeMap = permute(edgeMap, [3 1 2]);

% UNCOMMENT BELOW TWO LINES TO SAVE _detectionmaps.
% save([strcat(save_dir, '/', name) '_detectionmaps'], ...
%      'BWdilated', 'edgeMap');
 
% Provide notification if any raw_data rows contain all zeros. If there
% are zeros in one column, but not another, then do not generate a message,
% because the tracking function in MasterDataFrame.py should sort this.
[r, c] = find(raw_data(:, :, 1) == 0.0);
if numel(unique(r)) ~= numel(r)
    ['Saved .mat files for ', name, '. CHECK OUTPUT!']
else
    ['Saved .mat files for ', name]
end



    
    




% #########################################################################
% #################### Start of helper functions ##########################
% #########################################################################

function[min, max, xmin, xmax, ymin, ymax, letter, save_dir, comments] = ...
                                                    parse_inputs(varargin)
% Function to set the min and max frames, the region of interrest, 
% save letter, save directory, and comments saved based on arguments given.
% If the variable max = -1, the max will be set to the end of the video 
% after this function is called.
min = 1;
max = -1;
xmin = 1;
xmax = 512;
ymin = 1;
ymax = 512;
letter = '';
save_dir = pwd;
comments = '';
if(nargin == 1 && ~ischar(varargin{1}))
    min = varargin{1};
    max = varargin{1};
elseif(nargin == 1 && ischar(varargin{1}))
    comments = varargin{1};
elseif(nargin == 2 && ~ischar(varargin{2}))
    min = varargin{1};
    max = varargin{2};
elseif(nargin == 2 && ischar(varargin{2}))
    min = varargin{1};
    max = varargin{1};
    comments = varargin{2};
elseif(nargin == 3)
    min = varargin{1};
    max = varargin{2};
    comments = varargin{3};
elseif(nargin == 6)
    min = varargin{1};
    max = varargin{2};
    xmin = varargin{3};
    xmax = varargin{4};
    ymin = varargin{5};
    ymax = varargin{6};
elseif(nargin == 7)
    min = varargin{1};
    max = varargin{2};
    xmin = varargin{3};
    xmax = varargin{4};
    ymin = varargin{5};
    ymax = varargin{6};
    letter = varargin{7};
elseif(nargin == 8)
    min = varargin{1};
    max = varargin{2};
    xmin = varargin{3};
    xmax = varargin{4};
    ymin = varargin{5};
    ymax = varargin{6};
    letter = varargin{7};
    save_dir = varargin{8};
elseif(nargin == 9)
    min = varargin{1};
    max = varargin{2};
    xmin = varargin{3};
    xmax = varargin{4};
    ymin = varargin{5};
    ymax = varargin{6};
    letter = varargin{7};
    save_dir = varargin{8};
    comments = varargin{9};
end


function[image] = accumulate_frames(video, num_frames, place_in_series, num_accumulated)
% Function to accumulate adjacent frames to aid detection of weak signals.
% Adds together num_accumulated frames, centered on place_in_series. For 
% example, num_accumulated = 5 means that the program will accumulate 5 
% frames, 2 frames on each side of the frame number place_in_series. Note:
% num_accumulated must be odd! Also, flips image rightside-up since 
% fitsread flips it upside-down.

% Set min and max frame bounds:
ac_min = floor(place_in_series - (num_accumulated - 1) / 2);
ac_max = ceil(place_in_series + (num_accumulated - 1) / 2);

% Perform the addition:
if(ac_min < 1 && ac_max < num_frames + 1)
    image = sum(video(:, :, 1:ac_max), 3);
elseif(ac_min > 0 && ac_max > num_frames)
    image = sum(video(:, :, ac_min:num_frames), 3);
elseif(ac_min < 1 && ac_max > num_frames)
    image = sum(video(:, :, :), 3);
else
    image = sum(video(:, :, ac_min:ac_max), 3);    
end

% Flip image:
image = flipud(image);

% Get rid of super high intensity pixels:
background = initial_background_intensity(image);
temp_image = background_subtract(image, background);
rms = initial_background_rms(temp_image);
image(find(temp_image > 14 * rms)) = background;
% imshow(image, [min(min(image)) max(max(image))]);
% pause;


function[background] = initial_background_intensity(image)
% Function to define initial background intensity.
background = sum(sum(image)) / (size(image, 1) * size(image, 2));


function[image] = background_subtract(image, background)
% Function to subtract background after it has been calculated.
image = image - background;


function[smoothed] = gaussian_filter(image, gaussian_kernel_size, gaussian_std)
% Function to apply Gaussian filter of kernal size fspecial_kernel_size and
% standard deviation fspecial_std.
h = fspecial('gaussian', gaussian_kernel_size, gaussian_std);
smoothed = imfilter(image, h);


function[rms] = initial_background_rms(image)
% Function to calculate the RMS spread of the pixel values in an image.
meansq = sum(sum(image .^ 2)) / (size(image, 1) * size(image, 2));
rms = sqrt(meansq);


function[BWdilated] = clustering_detection(image, rms, threshold_RMS_multiplier, close_object, open_object, dilate_object)
% Function that detects clustering of high-intensity pixels. Performs 
% morphological processing after the detection.
level = threshold_RMS_multiplier * rms;
image(find(image < level)) = 0;
image(find(image ~= 0)) = 1;

BWc = imclose(image, close_object); %joins objects close together into a single object

BWc = imclose(BWc, strel('rectangle', [10 5])); %joins vertical pieces

BWco = imopen(BWc, open_object); %gets rid of dots of radius 2 or less
BWdilated = imdilate(BWco, dilate_object); %further dialates image (blurs it a bit)


function[background] = calculate_background_intensity(image, BWdilated)
% Refinement of background calculation based on first pass of clustering
% analysis.
ind2 = find(BWdilated == 0); % creates an index of all non-molecule points
background = sum(image(ind2)) / size(ind2, 1); %background is calculated using non-molecular points


function[rms] = calculate_background_rms(image, BWdilated)
% Refinement of background RMS calculation based on first pass of
% clustering analysis.
ind2 = find(BWdilated == 0);
meansq = sum(image(ind2) .^ 2) / size(ind2, 1);
rms = sqrt(meansq);


function[edges, edgeMap, connected, number_of_molecules] = edge_detection(image, BWdilated, edge_close_object, edge_open_object, edge_dilate_object)
% Function that performs the edge detection process. Perorms morphological
% processing after the detection.

% Identification of connected objects in picture. "number" is the number of
% molecules detected. "connected" is an array with 0 at background, 1s at 
% first molecule. 2s at second...
[connected, number_of_molecules] = bwlabel(BWdilated, 8);

% This option forces all the molecules to be classified as a single entity.
% number_of_molecules = 1;
% connected = BWdilated;

% Edge detection of molecular regions detected by clustering
edges = image;
[m, n] = find(connected == 0); %all pixels at which no molecule is detected by clustering
R = size(m, 1);
for s = 1:R; % cycles through each point using s as the index
    edges(m(s), n(s)) = 0;
end
edges = edge(edges, 'log');%finds the edges of molecules using the Laplacian of Gaussian method

edgeMap = imfill(edges, 'holes'); %fills regions enclosed by an edge, forming a map of molecules for edge detection
% imshow(edgeMap);
% pause;
edgeMap = imclose(edgeMap, edge_close_object); %same as before, joins close objects
% imshow(edgeMap);
% pause;
edgeMap = imopen(edgeMap, edge_open_object); %same as before, gets rid of noise (after the fill)
% imshow(edgeMap);
% pause;
edgeMap = double(edgeMap);
edgeMap2 = imdilate(edgeMap, edge_dilate_object); %same as before, further dilates the image (blurs it a bit)
% imshow(edgeMap);
% pause;
[m, n] = find(edgeMap2 == 0);
R = size(m, 1);
for s = 1:R;
    edges(m(s), n(s)) = 0; %Isolates molecular regions from empty spaces
end
% imshow(BWdilated);
% pause;
% imshow(edges * 100 + image, [min(min(edges * 100 + image)) max(max(edges * 100 + image))]);
% pause;


function[data] = catalog(image, connected, number_of_molecules)
% Cataloging of x- and y-coordinates and data compilation.
data = zeros(1, 6);
for i = 1:number_of_molecules;  %considers each molecule 1 at a time    
    [m, n] = find(connected == i); % m and n are the x-y values of pixels of molecule i
    R = size(m, 1);
    data_intensity = zeros(R, 1);
    data_x = zeros(R, 1);
    data_y = zeros(R, 1);
    for s = 1:R; % cycles through each point using s as the index
        data_intensity(s) = image(m(s), n(s));
        data_x(s) = n(s);
        data_y(s) = m(s);
    end
    
    % While still in the loop(that is considering each molecule at a time),
    %use the data to get info about each molecule.  Here is what it each
    %data point means:
    % data(i,1) = xCOM of molecule i
    % data(i,2) = yCOM of molecule i
    % data(i,3) = total intensity for molecule i
    % data(i,4) = Radius of gyration squared about x axis for molecule i
    % data(i,5) = Radius of gyration squared about y axis for molecule i
    % data(i,6) = overall Rg for molecule i
    % We use the approximation that intensity is proportional to mass.
    data(i, 3) = sum(data_intensity);
    data(i, 1) = sum(data_x .* data_intensity ) / data(i, 3);
    data(i, 2) = sum(data_y .* data_intensity ) / data(i, 3);
    data(i, 4) = sum(((data_x - data(i, 1)) .^ 2) .* data_intensity) / data(i, 3);
    data(i, 5) = sum(((data_y - data(i, 2)) .^ 2) .* data_intensity) / data(i, 3);
    data(i, 6) = sqrt(data(i, 4) + data(i, 5));
end


function[extension_lengths] = extension_length(edges, connected, number_of_molecules)
% Function to calculate extension lengths, given the edge detection. Simply
% finds the longest distance between any pair of points on the edge of a 
% molecule.
totalEdges = edges;
permConnected = connected;
extension_lengths = zeros(number_of_molecules, 1);
for i = 1:number_of_molecules
    connected = permConnected;
    maxLength = 0;
    isolatedEdges = totalEdges;
    %isolatedEdgeMap = edgeMap;
    connected(find(connected ~= i)) = 0;
    connected(find(connected ~= 0)) = 1;
    connected = imdilate(connected, strel('disk', 2));
    isolatedEdges(find(connected ~= 1)) = 0;%isolates a molecule for length calculation

    [m, n] = find(isolatedEdges ~= 0);
    R = size(m, 1);
    if(R > 0)
        tempCoordinates = zeros(2, 2);
        for s = 1:R;%find the longest distance between any 2 pixels on the edge
            for t = 1:R;
                tempLength = sqrt((m(s) - m(t))^2 + (n(s) - n(t))^2);
                if(maxLength < tempLength)%max distance between 2 points on edge
                    maxLength = tempLength;%max distance so far, will be length of molecule i
                    tempCoordinates(1, 1) = n(s);
                    tempCoordinates(1, 2) = m(s);
                    tempCoordinates(2, 1) = n(t);
                    tempCoordinates(2, 2) = m(t);
                end
            end
        end
        
%         xends = sort([tempCoordinates(1, 1) tempCoordinates(2, 1)]);
%         xs = xends(1):xends(2);
% %         display(size(xs));
%         ys = zeros(size(xs, 2));
%         for s = 1:size(xs, 2)
%             ys(s) = round((tempCoordinates(2, 2) - tempCoordinates(1, 2))/(tempCoordinates(2,1) - tempCoordinates(1,1)) * (xs(s) - tempCoordinates(1, 1)) + tempCoordinates(1, 2));
%             totalEdges(ys(s), xs(s)) = 1;
%             isolatedEdges(ys(s), xs(s)) = 1;
%         end
%         imshow(isolatedEdges)
%         pause

        extension_lengths(i) = maxLength;
    end
end


function[filtered_data] = filter_data(raw_data)
% Function to filter "molecules" that are too dim/too bright to consider.
%size(raw_data)
mean_ints = mean(reshape(raw_data(:, :, 3), 1, []));
%std_ints = std(reshape(raw_data(:, :, 3), 1, []));
min_thresh_int = mean_ints * 0.33;
max_thresh_int = mean_ints / 0.33;
filtered_data = raw_data;
for i = 1:size(filtered_data, 1)
    for j = 1:size(filtered_data, 2)
        if(filtered_data(i, j, 3) < min_thresh_int)
            filtered_data(i, j, :) = 0;
        elseif(filtered_data(i, j, 3) > max_thresh_int)
            filtered_data(i, j, :) = 0;
        end
    end
end
